{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For the 2-th left and right singular vectors and singular values, \n",
      "s * v = \n",
      "[-0.13491279  0.00735695  0.23851439  0.29752861]\n",
      "u * M = \n",
      "[-0.13491279  0.00735695  0.23851439  0.29752861]\n"
     ]
    }
   ],
   "source": [
    "# 验证第n组左右奇异向量与奇异值满足自洽方程\n",
    "\n",
    "dim = 4\n",
    "M = np.random.randn(dim, dim)\n",
    "U, S, V = np.linalg.svd(M)\n",
    "\n",
    "n = 2  # 检查第n个组左右奇异向量\n",
    "\n",
    "u = U[:, n]\n",
    "s = S[n]\n",
    "v = V[n, :]\n",
    "\n",
    "print('For the ' + str(n) + '-th left and right singular vectors and singular values, ')\n",
    "print('s * v = ')\n",
    "print(s * v)\n",
    "print('u * M = ')\n",
    "print(u.dot(M))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 最大奇异值及对应的左右奇异向量的迭代算法\n",
    "# 假设矩阵M为实矩阵\n",
    "\n",
    "def svd0(mat, it_time=100, tol=1e-15):\n",
    "    \"\"\"\n",
    "    :param mat: input matrix (assume to be real)\n",
    "    :param it_time: max iteration time\n",
    "    :param tol: tolerance of error\n",
    "    :return u: the dominant left singular vector\n",
    "    :return s: the dominant singular value\n",
    "    :return v: the dominant right singular vector\n",
    "    \"\"\"\n",
    "    dim0, dim1 = mat.shape\n",
    "    # 随机初始化奇异向量\n",
    "    u, v = np.random.randn(dim0, ), np.random.randn(dim1, )\n",
    "    # 归一化初始向量\n",
    "    u, v = u/np.linalg.norm(u), v/np.linalg.norm(v)\n",
    "    s = 1\n",
    "\n",
    "    for t in range(it_time):\n",
    "        # 更新v和s\n",
    "        v1 = u.dot(mat)\n",
    "        s1 = np.linalg.norm(v1)\n",
    "        v1 /= s1\n",
    "        # 更新u和s\n",
    "        u1 = mat.dot(v1)\n",
    "        s1 = np.linalg.norm(u1)\n",
    "        u1 /= s1\n",
    "        # 计算收敛程度\n",
    "        conv = np.linalg.norm(u - u1) / dim0 + np.linalg.norm(v - v1) / dim1\n",
    "        u, s, v = u1, s1, v1\n",
    "        # 判断是否跳出循环\n",
    "        if conv < tol:\n",
    "            break\n",
    "    return u, s, v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The dominant left singular vector by svd = \n",
      "[-0.20401265 -0.35018626  0.91419277]\n",
      "The dominant left singular vector by the iterative algorithm = \n",
      "[ 0.20401265  0.35018626 -0.91419277]\n",
      "\n",
      " The dominant right singular vector by svd = \n",
      "[ 0.03966432 -0.09903221  0.97139136 -0.00925571 -0.21197294]\n",
      "The dominant right singular vector by the iterative algorithm = \n",
      "[-0.03966432  0.09903221 -0.97139136  0.00925571  0.21197294]\n",
      "\n",
      " The dominant singular value by svd = \n",
      "2.51012165731717\n",
      "The dominant singular value by the iterative algorithm = \n",
      "2.5101216573171703\n"
     ]
    }
   ],
   "source": [
    "# 验证上面算法的正确性\n",
    "# 利用numpy库中的svd函数计算最大奇异值及对应的奇异向量\n",
    "M = np.random.randn(3, 5)\n",
    "U, S, V = np.linalg.svd(M)\n",
    "u0, s0, v0 = U[:, 0], S[0], V[0, :]\n",
    "\n",
    "# 利用上述算法计算最大奇异值及对应的奇异向量\n",
    "u1, s1, v1 = svd0(M)\n",
    "\n",
    "print('The dominant left singular vector by svd = ')\n",
    "print(u0)\n",
    "print('The dominant left singular vector by the iterative algorithm = ')\n",
    "print(u1)\n",
    "\n",
    "print('\\n The dominant right singular vector by svd = ')\n",
    "print(v0)\n",
    "print('The dominant right singular vector by the iterative algorithm = ')\n",
    "print(v1)\n",
    "\n",
    "print('\\n The dominant singular value by svd = ')\n",
    "print(s0)\n",
    "print('The dominant singular value by the iterative algorithm = ')\n",
    "print(s1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
